function [fitted_TAC] = fit_experimental_data(t_exp, A_exp, lambda, dt, k_0, its, T0, alpha, a, CA_ref)

%% Fit experimental time-activity data to the biokinetic model using SA
% t_exp, A_exp = experimental time-activity data for a given radiopharmaceutical 
% lambda = decay rate of the radiopharmaceutical
% dt = time step for TAC generation
% k_0 = initial set of parameters of the biokinetic model
% its = number of iterations to run (SA)
% T0 = initial temperature (SA)
% alpha = temperature modulation (SA)
% a, CA_ref = a*CA_ref weights the influence of the reference cumulated activity on the cost function

% Optimize k via SA
t_max = -log(0.01)/lambda;
t_array = 0:dt:t_max;
[k_opt, F_opt] = simulated_annealing(t_exp, A_exp,  k_0, CA_ref, lambda, T0, alpha, its, t_max, a, 0.1);

% Solve the model using optimal parameters
y_opt = model(k_opt, lambda, dt, t_max);

% Extract tumor time-activity data
A_opt =  y_opt(:,2)+y_opt(:,3);
A_opt(A_opt<0) = 0; 
fitted_TAC = [t_array; A_opt']';


end


function [k_best, F_best] = simulated_annealing(t_exp, y_exp, k_0, CA_ref, lambda, T0, alpha, its, t_max, a, dt)
    
% Find the optimal set of parameters k by minimizing F_cost using a
% simulated annealing (SA) algorithm


% array de tempos para F_custo
t_array = 0:dt:t_max;
for i = 1:length(t_exp)
    [~, aux] = min(abs(t_array - t_exp(i)));
    ind_array(i)=aux(1);
end

% initial solution and cost
y_best =  model(k_0, lambda, dt, t_max);
[F_best] = F_cost(y_exp, t_array, (y_best(:,2)+y_best(:,3)), CA_ref, ind_array, dt,a);
k_best = k_0;

% run simulated annealing
F_min = F_best;
k_min = k_0;
T = T0;
for i = 1:its
    for j = 1:its

        % randomly perturbate k
        k_i = k_best;
        ind = fix(rand * length(k_0)) + 1;
        k_i(ind) = k_i(ind) * (1 + (rand - 0.5) * T / T0);
       
        % solve model for the new k and compute cost
        y_i = model(k_i, lambda, dt, t_max);
        [F_i] = F_cost(y_exp, t_array, (y_i(:,2) + y_i(:,3)), CA_ref, ind_array, dt, a);
       
        % accept new solution
        if F_i < F_min || rand() < exp(-(F_i - F_min) / T)
            F_min = F_i;
            k_min = k_i;
            
            % update best solution
            if F_min < F_best
                F_best = F_min
                k_best = k_i;
            end

        end
    end
    T = alpha * T;
end
end

function [y] = model(k, lambda, dt, t_max)

%     Biokinetic model.
%     Compartments (5): 
%     y1 = Blood
%     y2 = Tumor 1 
%     y3 = Tumor 2 
%     y4 = Remaining Tissues 1
%     y5 = Remaining Tissues 2
%     Coupling constants (7):
%     k1 = Blood -> Tumor 1 
%     k2 = Tumor 1  -> Blood
%     k3 = Tumor 1  -> Tumor 2
%     k4 = Tumor 2 -> Tumor 1
%     k5 = Blood -> Remaining Tissues 1
%     k6 = Remaining Tissues 1 -> Blood
%     k7 = Remaining Tissues 1 -> Remaining Tissues 2


% Solve model using Euler's method
t=0:dt:t_max;
y(5, length(t))=0;
y(1,1)=1;
for i=2:length(t)

    y1 = k(6)*y(4, i-1) + k(2)*y(2, i-1) - ( k(1) + k(5) + lambda )*y(1, i-1);
    y2 = k(1)*y(1, i-1) + k(4)*y(3, i-1) - ( k(2) + k(3) + lambda )*y(2, i-1);
    y3 = k(3)*y(2, i-1) - ( k(4) + lambda )*y(3, i-1);
    y4 = k(5)*y(1, i-1)- ( k(6) + k(7) + lambda )*y(4, i-1);
    y5 = k(7)*y(4, i-1) - lambda*y(5, i-1);
    
    y(1,i) = y(1,i-1) + y1*dt;
    y(2,i) = y(2,i-1) + y2*dt;
    y(3,i) = y(3,i-1) + y3*dt;
    y(4,i) = y(4,i-1) + y4*dt;
    y(5,i) = y(5,i-1) + y5*dt;

    ind=y<0;
    y(ind)=0;
end
y = y';

end



function [cost] = F_cost(y_exp, x_teo, y_teo,  CA_ref, ind_array, dt, a)
        
         % Cost function. For a~=0, cost is modified using the reference
         % cumulated activity CA_ref (weighted using a)
         CA_teo = trapz(x_teo, y_teo); 
         ind=fix(2/dt);
         for i=1:length(y_exp)
             ind2=max(1, ind_array(i)-ind):(ind_array(i)+ind);
             aux3=abs(y_teo(ind2)-y_exp(i))'; 
             aux2=abs(ind2-ind_array(i))*dt;
             aux(i) = min( aux3.^2 / (max(y_exp)^2)  + 0.001*aux2.^2);
         end
         s1 = sum(aux);
         s2 =  a*((CA_teo-CA_ref)^2/CA_teo^2);
         cost = s1 + s2;
end